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Abstract 

The algorithm’s objective is to efficiently solve Dynamic Linear Programs (DLP) by taking advantage of their 
special staircase structure. This algorithm constitutes a stepping stone to an improved algorithm for solving Dynamic 
Quadratic Programs, which, in turn, would make the nonlinear programming method of Successive Quadratic 
Programs more practical for solving trajectory optimization problems. The ultimate goal is to bring trajectory 
optimization solution speeds into the realm of real-time control. 

The algorithm exploits the staircase nature of the large constraint matrix of the equality-constrained DLPs 
encountered when solving inequality-constrained DLPs by an active set approach. A numerically-stable, staircase QL 
factorization of the staircase constraint matrix is carried out starting from its last rows and columns. The resulting 
recursion is like the time- varying Riccati equation from multi-stage LQR theory. The resulting factorization increases 
the efficiency of all of the typical LP solution operations over that of a dense matrix LP code. At the same time 
numerical stability is ensured. The algorithm also takes advantage of dynamic programming ideas about the cost-to-go 
by relaxing active pseudo constraints in a backwards sweeping process. This further decreases the cost per update of 
the LP rank-1 updating procedure, although it may result in more changes of the active set than if pseudo constraints 
were relaxed in a non-stagewise fashion. The usual stability of "closed-loop" Linear/Quadratic optimally-controlled 
systems, if it carries over to strictly linear cost functions, implies that the savings due to reduced factor update effort 
may outweigh the cost of an increased number of updates. 

An aerospace example is presented in which a ground-to-ground rocket's distance is maximized. This example 
demonstrates the applicability of this class of algorithms to aerospace guidance. It also sheds light on the efficacy of 
the proposed pseudo constraint relaxation scheme. 


Introduction 


The objective of the present work is to develop and test a special-purpose algorithm for the solution of 
Dynamic Linear Programs (DLP). The general form of a DLP is as follows: 
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where the x, vectors constitute the decision vector time history, the c, vectors are linear cost coefficients, and the A Li 
and Au +1 matrix blocks and the b, vectors define the linear problem constraints. The bracketed equality and inequality 
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signs in eq. lc indicate that both forms of constraints may be present; some rows may be equalities while others arc 
inequalities. The problem in eq. la-lc is also known as a staircase LP. 

A reason for interest in this problem from an aerospace controls point of view comes from its relationship to 
the following multi-stage trajectory optimization problem: 


find: 
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to minimize: 

N-l 

J = £ L(x k ,u k Jc) + <(»(x N ) 

k=0 


(2b) 

subject to: 

x 0 given 


(2c) 


x k+i = f ( x k» u k* k ) 

fork = 0 ... N-l 

(2d) 


a(x k ,u k dc) 0 

for k = 0 ... N-l 

(2e) 


a(x N ,N) {<} 0 
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where the Uj and x i+1 vectors constitute the control and state vector time histories, the LO and <|>0 functions define the 
nonlinear stagewise and terminal cost s, eq. 2c defines the initial conditions, eq. 2d is the discrete-time dynamics 
difference equation, and constraints such as 2e and 2f may be present to restrict the states, or the control inputs, or 
both. 


The current paper is part of a research program that seeks a fast and reliable way to solve the problem in eq. 2a- 
2f. The ultimate goal is to do real-time aerospace guidance by repeatedly solving this problem. The program is 
taking a two-pronged approach, algorithm improvement and parallelization of computations. This paper relates to the 
first prong, algorithm improvement. Fletcher’s penalty function, trust region adaptation of the method of 
successive quadratic programs (SQP) [1] is one algorithm for solving such a nonlinear program (NP). This algorithm 
has fast local convergence properties, it ensures global convergence (to a local minimum), and it is good at handling 
inequality constraints. The application of this algorithm to the nonlinear trajectory optimization problem results in 
Quadratic Programming (QP) sub-problems with a special structure, the dynamic programming structure. Efficient 
solution of the NP requires efficient solution of the dynamic QP (DQP). Special-purpose algorithms for efficient 
solution of a DLP can be similar to special-purpose algorithms for efficient solution of a DQP. Thus, the present 
paper, in concentrating on DLP, constitutes a sort of warm-up exercise for later development of a DQP algorithm. 

In addition to providing a warm-up, the present work provides a point of comparison with other research efforts 
in the field. Little or no work has been done on special-purpose DQP algorithms [2,3], but much attention and effort 
has been devoted to special-purpose DLP algorithms (e.g., Refs. 4-7). Fourer provides a useful overview of different 
avenues of approach that have been tried [7]. He groups algorithms for problem la-lc into three categories: Compact 
Basis, Nested Decomposition, and Transformation. All get the correct answer, but none have proved particularly 
successful in that none consistently out-perform the general sparse simplex method with regard to computation time. 

The present algorithm is in the compact basis category; it works with a staircase factorization of the active 
constraints. The factorization used, a staircase QL factorization, is consistent with the plan for subsequent upgrading 
to handle the quadratic cost case. It has numerical stability, and there is no trade-off between numerical stability and 
factor compactness; general sparse matrix LP codes must deal with such trade-offs. The focus of the entire project is 
on aerospace guidance problems, hence the submatrices of the problem, the and A^blocks, are relatively dense. 
Therefore, there is hope that the current algorithm will out-perform general sparse matrix LP algorithms on these 
problems (e.g., algorithms such as MINOS [8]). 

The body of this paper concentrates on explanation of the algorithm, with an example and conclusions at the 
end. Before describing the algorithm, problem la-lc is related to a general LP on the one hand and to a control-type 
LP on the other hand. The algorithm description begins with a review of the application of the L! penalty function 
method to a general LP. The staircase QL factorization then gets presented along with methods for multiplier 
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solution, decision vector solution, and rank-1 update. The algorithm explanation concludes with the presentation of a 
specialized order for problem solution that could further reduce the computational burden. The example at the end of 
the paper demonstrates the algorithm's applicability to aerospace trajectory optimization and examines whether the 
special solution ordering yields increased computational efficiency. 

Equivalent Problem Forms 

The problem in eq. la-lc is a special case of the following general LP form: 


find: 

X 

(3a) 

to minimize: 

J = C T X 

(3b) 

subject to: 
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where x is defined in 

eq. la and A, b, and c are defined as: 
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Form 3a-3c is fully general. It is the LP form used in developing general active constraint algorithms [1]. 
A more specialized DLP problem statement clarifies the relationship of these problems to controls: 


find: 
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x k+ i = F k x k + G k u k + h k 


for k = 0 ... N-l 

(5d) 


A x *k + A u “k-b x » 1 

k k k ' 

r-i 

L<J 

[ 0 fork = 0... N-l 

(5e) 


A x N *N b XUN { = } 0 



(5f) 


where there is a direct correspondence between eq. 2a-2f and eq. 5a-5f, all functions having only linear and constant 
terms in the latter problem. The following definitions put problem 5a-5f in the format of problem la-lc: 
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Thus, the problem in eq. la-lc is related to controls. 


for k = 1 ... N-l 
for k = 0 ... N-l 
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Algorithm 


The Lj Exact Penalty Function and an Active Set LP Method* 

The L! exact penalty function method enforces equality and inequality constraints as in 3c by adding a penalty 
cost to the problem cost that is a weighted L, norm of the constraint violations. The penalty function reformulation 
of problem 3a-3c becomes 

find: X (^ a ) 

to minimize: J = c T x + x - b e ll, + lltA^ x - bjnl^llj } (7b) 

where |imu is a large positive constant, where constraints [A,b] in eq. 3c have been split into the equality constraints, 
[A e ,bJ, and the inequality constraints, [A m ,b m ], and where [a] + = max(a,0). This technique is actually very similar to 
the adjoining of constraints in Lagrange and Kuhn-Tucker formulations. The only difference is that the constraint 
multipliers are effectively limited in magnitude by tw If IW is greater than the magnitude of the largest multiplier 
in the solution of problem 3a-3c, then problem 7a-7b has the same solution. The advantage of this technique is that it 
solves the problem in a single phase, optimizing while achieving feasibility. It is a variant of the LP technique 
known as the big M method. The primary reason for using it in the current paper is to make the algorithm 
compatible with the proposed method for eventually solving the NP in eq. 2a-2f. 

Active constraints refer to those constraints that are satisfied as strict equalities. The active set method for 
solving problem 7a,7b consists of the following steps. 


1. Guess solution, x, and split [A.b] into active and inactive rows: 


= P [A.b] 

-^in*ct*^in»ct- 


where P is a permutation matrix, and where the active constraints at the guessed solution must 
yield an A^ that is square and nonsingular. 

2. Compute A**' 1 

3. Compute u.* = - (A.L) T (c + A T iMC 0i inact ) where the elements of lkn.ct are either 0, or +1W 


depending on whether the corresponding inactive constraint is an equality or an inequality and 
depending on whether it is positive or negative at the guessed solution, x. 

4. Find the element, i, of u*., that is furthest out of the allowable range for the penalty function, 

(U«t)i e Hw, +p nu ,J for equality constraints, (u«*)i e [0, +1W] for inequality constraints. Stop 
if no elements are out of range; the current x is optimum. 

5. Compute the search direction, Sx = A,*^, where ei is the unit vector with all zeros except for a 1 


in row i. 

6. Compute the new guessed solution, x = x + a Sx, where a has the same sign as (u lct )j and where 
its magnitude is chosen to be the smallest value that makes one of the inactive constraints active, 
min loti such that [A in , ct (x+ a Sx) - b^Jj = 0, for some row j. 


t The discussion of this section deals with the problem form in eq. 3a-3c. The discussion is based 
on algorithms presented in Ref. 1. 
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7. Interchange rows i and j between [ A^.b^J and tA inict ,b iniel ], update A, cl , and go to step 3. 


The equation in step 3 results from differentiation with respect to x of the augmented cost in eq. 7b. Despite the 
nondifferentiability of eq. 7b for some values of X, differentiation can be done if the active constraints are treated as 
adjoined equality constraints with unknown multipliers U. CI rather than as L] penalty terms. Steps 4 and 5 determine a 
descent direction. The step length is determined in step 6. The step length rule ensures that the entire step is in a 
descent direction of the piecewise-linear augmented cost. Because of the step length rule, the new active set changes 
by only one row. Step 7 takes advantage of this fact in the recomputation of A tct 's inverse. 

The algorithm used in this paper starts out by assuming that a set of pseudo constraints are active. This yields 
the identity matrix for the initial A, ct , and step 2 is trivial. The allowable range for the pseudo constraint multipliers 

is different than for the actual problem constraints, e [0,0]. They get dropped from the active set in the course 
of the algorithm unless the problem solution is not unique. 

For each constraint addition/deletion the algorithm cycles through steps 3-7. Most of the computational load 
per cycle is caused by manipulations with the inverse of the A iCt matrix, multiplication by it in steps 3 and 5 and 
rank-1 updating of it in step 7. The main idea of this paper — indeed the main idea of all compact basis schemes - is 
to compactly represent factors of A tcl that are suitable for carrying out multiplications by A tc /s inverse and that are 

easy to update when A act undergoes a rank-1 row change. 

Staircase QL Factorization for Staircase LP 

A acl inherits a staircase structure from A as in eq. 4. The compact QL factorization used here performs a stage- 
wise backwards sweep to factor the nonzero blocks of A act . The following recursion yields matrices that constitute a 
staircase QL factorization of A act : 


Dnn = Ann 

act 

Dk-lk-1 0 
Dkki L 

QooDoo = Lqo 


Ak-lk-l A k _ lk 

a. i*. i <ct * 

Qkk 0 D kk 


for k = N, 


.1 


(9a) 

(9b) 

(9c) 


where the A k . lk . llict and A k . lktct matrices are the nonzero blocks of the staircase A act matrix, and where the are 
orthonormal and the are lower triangular. The QL factorization is stored in the matrices Qj^, L^, and for k = 
N.....0 and the matrices D^., for k = N.....1. The matrices have no special properties. Equations 9b and 9c are 
only implicit relations for these factors, but the factors can be explicitly evaluated via Householder transformations. 
At stage k, the factorization begins with the data A k _i k _ liCt , A k . lktct , and D^, and it computes Q^, L^, D k .j k .i, and 

D^ !. The result D k . lk .! then completes the necessary data for stage k-1. 

Numerical stability of the factorization is ensured by the use of orthogonal transformations only. The 
computational complexity of the factorization algorithm is linear in the number of stages and cubic in the 
dimension(s) of the blocks, which is as efficient as can be expected if the blocks are dense. The factor storage is linear 
in the number of stages and quadratic in the dimension(s) of the blocks, which again is the best that can be achieved 
with dense blocks. 

Equations 9a-9b are a DLP equivalent to the matrix Riccati equation of time-varying, multi-stage Linear 
Quadratic Regulator theory. The lower blocks of the right hand side of eq. 9b act as a closed-loop dynamic difference 
equation as will be shown in the next section. The upper block on the right hand side, [D k . lk _i, 0], propagates active 

constraint effects backwards in time; it summarizes the constraints that x*.i must satisfy in order to make possible the 
satisfaction of all constraints from stage k-1 onwards. Note that the number of rows in [D k . lk .!, 0] does not 
necessarily equal the number of rows in [A k _i k _ licl , A k . lkact ]. 
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Staircase QL Solution for the Multiplier and Decision Vector Time Histories 

The basic operations involved to do steps 3 and 5 of the LP algorithm presented above involve solution of a 
linear system by orthogonal transformation and forward or backward substitution. This is similar to the forward and 
backward substitutions of the simplex method and its variants. The transformations and substitutions are done in 
stage-wise recursions using the stagewise factors. Performing the following backward recursion then forward recursion 
yields the active constraint multipliers. 

Backward recursion for intermediate multipliers ^ through 2^ 

l nn hi = - c N - A »w<*w ' (10a) 

K = - I>Lk K+I - c k - An.rfU'w ■ for k = N- 1 , .... 1 (10b) 

k>00 lo = -Du>h " C 0 " (Ifc) 


Forward recursion for intermediate multipliers Bo through and for active constraint multipliers through u N : 

act act 


) = Qooio 


(Ha) 


_ fik+1 


T 

: Qk+lk+1 


L \+i 


for k = 0 ... N-l 


(Hb) 


Un„=6, 


(lie) 


Step 5 of the LP algorithm is accomplished by solving the system of equations A ac ,5x = e,. To illustrate how 
the staircase QL factorization does this, the following equations present its use in solving the alternate system A acl x 
= b act . Again, a backward recursion followed by a forward recursion yields the solution. 


Backward recursion for intermediate nonhomogeneous constraint terms d N through d 0 and g N through g 0 : 

(12a) 




A. I 

^k-1 

= Qkk 

* »ct 

gk 



for k = N 1 


So “ Qoo do 

Forward recursion for the decision vectors x 0 through x^,: 
L 0 x 0 = go 

Lki x k = - D^x^, + g k 


for k = 1 ... N 


(12b) 

(12c) 

(13a) 

(13b) 


As stated earlier, eq. 13b is like the closed-loop dynamic difference equation of multi-stage LQR theory. The matrix 
Lu is lower triangular and allows for easy solution for Xj. in terms of x k ] and g k . 


Rank-1 Update of Staircase QL Factorization 

This section explains how to efficiently update the staircase QL factorization after a single constraint 
addition/deletion. This procedure must be carried out every time step 7 of the main algorithm is encountered. One 
could recompute the entire factorization, but the practicality of all LP codes hinges on their ability to update the 
factors for much less work than would be required to recompute them from scratch. 
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The general add/drop updating scheme for the staircase QL factorization must update the results of eq. 9a-9b 
when an arbitrary row j at stage kadd gets added to the active constraint set and another arbitrary row i at stage kdrop 
gets deleted from the active constraint set. Thus, [A^, A kk+liCt ]l k « kadd gets a new row and [A kkact , Akk+ 1 ict l tk=kdrop 
loses a row. The stages kadd and kdrop can have any relationship to each other, and the update algorithm must be able 
to handle all possible cases. Three different cases can occur, kadd > kdrop, kadd = kdrop, and kadd < kdrop. 

Efficient rank-1 update can be accomplished by a series of stage-wise rank-1 updates linked together in an 
appropriate manner. Three different stagewise rank-1 updating algorithms are needed to do this. The first algorithm 
updates the factors computed in eq. 9b when a new row has been added to the bracketed expression on the left-hand side 
of that equation. The second algorithm updates these same factors in the case of a row deletion from the bracketed 
expression on the left-hand side. This second algorithm also modifies Q^.j by a single Householder transformation. 
The third stagewise algorithm updates these same factors when has undergone an arbitrary rank-1 change. Recall 
that the bracketed matrix on the left-hand side of eq. 9b represents the input data for a given stage and the matrix 
together with the bracketed expression on the right-hand side represents the result of the stagewise factorization. The 
following discussion explains each of these three algorithms and the way in which they work together to accomplish 
the multi-stage rank-1 update. 

First, consider what happens to the stage k factorization when a new row gets added to either [A k _ lk _ ltcl , A k _ lkacl ] 
or Dj^. The algorithm begins by adding a row and a column to with all Os except for a 1 at the intersection of the 
new row and the new column. Thus, remains orthonormal. Suppose the new constraint row is [p T kl , p T k2 ], then 
the new row and column of are added so that eq. 9b temporarily becomes: 


Qum 

0 

Qkki2 

At. i k . i 

‘act 

Ak - lk .ct 


Dk-lk-l 

0 

o T 

1 

0 T 

P T kl 

P T k 2 

= 

p T kl 

P T k 2 

l 

e 

N> 

0 

Qkk 22 

0 

Dkk _ 


_ D kk -i 

k*kk _ 


where the Q kki j matrix blocks are just the blocks of the original matrix. Suppose n k is the dimension of the % k 
decision vector. Then it is also the dimension of the square lower-triangular matrix L^. A series of n k Givens 
rotations can be performed to zero out p T k2 while preserving the lower-triangular structure of L^. The first Givens 
rotation uses the last row of as the pivot row and zeros out the last element of p T k2 , and successive rotations use 
successively higher rows of L kk as the pivot and zero out successive elements of p T k2 going from right to left. 
Suppose these rotations are Gi to G nk . Then the new stage k factorization becomes: 


Qkknew - G "k‘ 


•Gi 


Qkkn 0 

• 0 T 1 

_ Qkk21 0 


Qkki2 

0 T 

QkJc22 . 


(15a) 


Dk-lk-t 

0 - 


D k -i k -i 

o * 

cPk, 

0 T 

= 

G nk *...*Gi* 

P T ki 

P T k 2 

_ f^kk*l ncw 

V - 


. Dkk 1 

Lkk . 


r D k.ik 

-l " 



D lc-lk-l ncw 

ii 

i — 





(15b) 


(15c) 


where the last equation has been included to emphasize the fact that the new D k . lk .j differs from the old D k . lk . 1 by only 
a single new row. This fact sets the stage for the use of this same algorithm at stage k-1. The new Q kk is 
orthonormal because the augmented matrix is orthonormal and because all Givens rotations are orthonormal. Thus, 
the new factors have all of the required properties for use in the LP algorithm described above. 
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Next, consider what happens to the stage k factorization when a row gets deleted from either [A k _ lk _ 1#ct , A k _ lkact ] 
or D^. The following development is based on ideas for QP from Ref. 9. Write Q& in the form 


Qkk “ 


" Qkkil 

<J T kk21 


0kki2 Qkkl3 
T 

Okie 22 Q kk23 


L Qkk3i 0kk32 0^33 


(16) 


where the middle column conforms in matrix multiplication with the constraint row that is getting deleted -- rows of 
can be referred to as constraints; they are propagated active constraints. The bottom blocks, [Qkk 3 i> ( lkk 32 ’Qkk 33 l’ 

have n k rows, the same as in the bottom blocks on the right hand side of eq. 9b, 


The stagewise deletion algorithm starts with a Householder transformation in which the g^ 12 row in the above 
representation is used as the pivot row to zero out qtki 2 ^ rows. Next, a series of n k Givens rotations is used 

to zero out successive elements of q m2 starting with the topmost element and working downwards. Again, the 
row in the above representation is used as the pivot. If the Householder transformation is H and the Givens rotations 
are G x to G nk , then the following changes to the stage k factorization result: 


Qwm.new ^ Okkl2new 
0 T 1 0 T 
. Okk21ncw ® Qttc22ncw 


= G nk ...G 1 -H 


Qkkil *lkkl2 Qkk 13 
q^kk21 ^btk22 Q T kk23 
L Okk3i Qkk32 Okk33 


Qkkncw 


^^llnew Q^^ncw 
0^21 new Q^22new 


^-lk-l n< 

P T kl 


■*kk- 


0 

P T k2 


rDk-ik-i o i 

= g i4 -..,g 1 .h{ Dkk _ t Lkk J 




(17a) 


(17b) 


(17c) 


where [p T kl ,p T k2 l corresponds to the constraint that is getting dropped. Orthonormality of the original matrix 
ensures the form of the result on the left-hand side of eq. 17a. Note that the matrix D k _ lk .i oew is a function only of 
D k . lk _i and H; the Givens rotations do not affect it. Therefore, another Householder transformation, H", can be 
constructed based on the same Householder vector. It yields: 


Dk 

dV 


lk-ln 


1 drop 


H" D k . lk _i 


(18) 


where d k . ldrop is not necessarily equal to p kl . This sets the stage for propagation of the constraint deletion process 
backwards to stage k-1. If Qk.i k .i gets transformed according to 


Q k 


Ik- 1 interim 



(19) 


then Qk-ik-iintenn, is still orthonormal because H" is orthonormal, and because H” is equal to its transpose, constraint 
[0 T ,d T k _ ldrop ] is the constraint that must get dropped at stage k-1. The foregoing algorithm can accomplish the deletion 
at this next preceding stage. 
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The algorithm that performs the multi-stage rank-1 update of the staircase QL factorization starts with the 
highest stage at which either a constraint addition or deletion occurs. It uses whichever of the two foregoing stagewise 
updating algorithms is appropriate to propagate the addition or deletion backwards. It continues until it reaches a stage 
at which both an addition and a deletion must take place. One, but not both, of the changes at this stage may be the 
result of a backwards propagation. At this stage of the concurrent add/drop, the multi-stage algorithm first does a 
single-stage constraint addition followed by a single-stage constraint deletion with no change of stage in between. 

If the index of this stage is k, then D k . lk .j will differ from its pre-update value by a rank-1 change at most. 
This can be shown by recognizing that the result on the left-hand side of eq. 15c becomes the input data on the right- 
hand side of eq. 18 when an add followed by a drop both occur at the same stage: 


Du _l k i 
* l ncw 

_ dT k-l dro p 


= H" 


D k lk-1 " 
. ^ki _ 


( 20 ) 


H" is a Householder transformation; it differs from the identity matrix by a rank-1 matrix, hence the conclusion about 
the change in Define this rank-1 change in terms of the vectors r k _ t and s^: 

Dk-lk-l ncw = D k .i k _! + r^sVl (21) 


If either r k . x or s k _j is the 0 vector, then the multi-stage rank-1 update is complete. If not, then another stagewise 
updating algorithm is needed. 


The final stagewise updating algorithm must update the stagewise factors for an arbitrary rank-1 change in the 
data D^. It is allowed to produce at most a rank-1 change in D klkl . This restriction on its effect on D k . lk ., makes it 
self recursive for all subsequent stagewise factorizations in the backwards chain. It can be used for updating the 
factorizations of all stages that precede the concurrent constraint addition/deletion stage. It can be used recursively 
until no more updating is needed. 


One might suppose that the necessary algorithm has already been developed in a work such as Ref. 10. That 
paper is a good reference for rank-1 modifications, and it defines the general methodology used in the algorithm below, 
but the relevant algorithm from [10] would result in a rank-2 change to D k _ lk _!. This would destroy the stagewise 
recursive applicability of the algorithm, hence the modified algorithm presented below. 


Suppose there has been a rank-1 modification to as in eq. 21 (except at stage k instead of stage k-1). Then, 
eq. 9b gets modified: 


Qkkj 

where 


A k -ik-i lcl A k . lkjct 

0 [^kk^'l'kS^'k] 


Dk-lk-1 0 
Dkk-1 4k 



( 22 ) 


* Vk -* 1 ~ 

- Qkk 

■ 0 ■ 

L w k J 

_ r k _ 


(23) 


and where v k _j and D k _ lk _i have the same number of rows, n dk] . The algorithm starts by reducing the v-w vector to a 
vector with zeros in all of its entries except the last two. This is done by first applying a Householder transformation, 
H lt to the first n dkl +l rows to zero out the first n dM rows. Then a series of Givens rotations, G, to G^, is applied 

to successive pairs of rows of the resulting vector to zero out successive elements until only the last two elements are 
left nonzero. These same transformations are applied to all terms on both sides of eq. 22, and the two terms on the 
right hand side of the equation are added together with the following (partial) result: 
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(24) 


v k _, 0 0 ... 0 “ 

* * 

* * * 0 

* * * * 

* * * * * 

* * * * * — 

where a is a scalar and where asterisks (*) indicate nonzero scalar elements of the matrix. The top row of vector 
entries in the right hand matrix corresponds to the upper right-hand 0 block in the bracketed expression on the right 
side of eq. 9b. The bottom rows constitute a matrix with nonzero elements on the first diagonal above the main, on 
the main diagonal, and below the main diagonal. All elements on diagonals that are 2 or more above the main 
diagonal are zero. The nonzero entry in the first column of the top block, av k .j, results from application of the H ; 
Householder transformation to matrix [O^L 1 ^] 7 . 

The remaining transformations are applied in order to restore lower triangularity to the matrix on the right hand 
side of eq. 24. First, a series of n k -l Givens rotations, G nk to G^.* is applied to successive pairs of rows of the 

matrix starting from the last two rows and working up to the first two rows in the lower block. Each rotation zeros 
out one of the above-diagonal elements. At the end of this operation the matrix has the form 




l* T k - 


on 


otv k _i 0 0 ... 0 “ 

* 

* * 0 
* * * 


(25) 


* * * * I 

* * * * * 


so that the lower block is lower triangular. The final part of the algorithm is to apply a last Householder 
transformation, H 2 , to zero out the first column of the top block. These operations result in the following factor 

updates: 

Qkw = H 2 *G 2nk . 2 *... , G 1 *H 1 *Q kk (26a) 


|_ ^^new ^^new 

In both of the Householder transformations, the first n^.j elements of the transformation vector are parallel to v k . b and 
none of the Givens rotations affect the first n dk _, rows of the bracketed expression on the right of eq. 26b. Therefore, 
D k -ik-i ncw differs from D k _ lk .! only by a rank-1 matrix: 


= H 2 *G 2nk - 2 v..-G 1 .H 1 | ^ L J + [ “][ ° T s T k ]} (2®) 


Dk-ik-i^ = Dk-ik-i + v k -iy T k-i 

where the vector y k .j can be determined from the algorithm presented above. Thus, the algorithm updates stage k 
according to the rank-1 change in the stage's input data, and it produces a similar rank-1 change in the input data for 
stage k-1 The multi-stage rank-1 updating algorithm propagates these rank-1 changes backwards until at some stage 
one or both of the vectors in the rank-1 change are zero. This occurs at least by the time stage k - 0 is reached 
because D.i.j has zero dimension. 
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Numerical stability of the factorization update is ensured by the use of orthogonal transformations only. The 
computational complexity of the multi-stage update algorithm is linear in the number of stages affected and quadratic 
in the dimension(s) of the blocks, which is as efficient as can be expected if the blocks are dense. If the number of 
stages affected by a particular row interchange of active and inactive constraints can kept small, then the cost of the 
update will be small. This fact provides the motivation for the solution scheme presented in a later section. 

Possible Improvements to Banded Staircase QL Factorization 

Several issues come to mind in considering the forgoing use of a QL factorization for an LP basis factorization. 
They all revolve around a single question: is the entire factorization needed to implement the LP algorithm? For 
instance, a general LP method has been developed that uses LQ factorization but does not store Q [11]. Not storing 
the Q factors would yield a great savings in memory and computation time if it carried over to the present multi-stage 
algorithm. This presents no difficulty to the procedures for solving for the multiplier and decision vectors, steps 3 and 
5 of the main LP algorithm. The problem with not storing Q occurs in the factor update, step 7. There is no 
apparent way to do the single-stage constraint deletion or the single-stage rank-1 modification without storing at least 
some of the matrix. Reference 9 has some ideas in its section on quadratic programming that could be used to 
eliminate storage of the lower part of Q^. Alternatively, storage of and could be eliminated. Savings in 

computation time and memory would be about the same for either scheme, about 30% savings. These issues may be 
explored in a later work. 

Backwards-Sweeping Pseudo Constraint Relaxation and an Alternate Method of Selecting the 
Active Constraint to Drop 

In theory, all dynamic programming problems can be solved by first computing the cost-to-go at each stage, 
then solving a single stage optimization at each stage. Part of the cost for each of these single stage problems is the 
cost-to-go that results from the stage’s decisions. For DLPs and for their associated Lj penalty function problems, the 
cost-to-go at a given stage is a piecewise-linear convex function of the decisions at that stage. This convexity 
property gives rise to a hope that DLPs may have a property like the stability property of their quadratic-cost 
counterparts, multi-stage LQR problems. In the DLP context, this property might mean that a small change in the 
decisions at a given stage would give rise to even smaller changes in the state at subsequent stages. This might 
translate into a grouping of constraint additions and deletions at stages nearly following the stage at which the decision 
variations are taking place. 

If this property exists, it can be exploited without the necessity of computing the entire cost-to-go function. If 
all of the active constraint multipliers for constraints following a given stage are within their allowable range, then the 
guessed solution is an optimal trajectory for all stages following that stage. Also, the local linear piece of the cost-to- 
go function is known. Suppose the given stage can be optimized without causing any of the multipliers at subsequent 
stages to exceed their penalty function bounds. Suppose also that all of the original pseudo constraints are active 
for the preceding stages. Then, the rank-1 updates that would have to be done during the optimization of that stage 
might involve changes to very few stages. The assumption about the the pseudo constraints ensures that the updates 
will not affect any of the stages preceding stage k-1 if stage k is being optimized. The possibility of stability implies 
that max(kadd,kdrop) might, in most cases, not be much larger than k. 

A change is needed to the main LP procedure presented above. It allows the multipliers at stages subsequent to 
the stage being optimized to vary outside of their penalty function bounds. The modification needs to be in the 
selection of the active constraint that gets dropped on each cycle. In the main algorithm, the dropped constraint is the 
same as the non-optimal constraint that gets relaxed in steps 4-6. This could cause an active constraint multiplier that 
was within its bounds to go out of its bounds. If the multiplier corresponded to a constraint at a subsequent stage, 
then the optimality of the subsequent stages would break down. 

This situation can be avoided by performing a search in the active constraint multiplier space for the active 
constraint to be dropped. This search is the dual of that carried out in steps 5 and 6 of the main algorithm, and the 
search direction is defined by relaxing the penalty function constraint on the multiplier associated with inactive 
constraint j, the inactive constraint that is becoming active. The size of the step in multiplier space is chosen to be 
the smallest that brings one of the active constraint multipliers to a bound which the multiplier would violate if the 
step size were larger; the new active constraint must be included in this test. The active constraint whose multiplier 
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bound limits this step size is the active constraint that gets dropped. It is not necessarily the constraint whose 
relaxation was dictated in step 3 of the main algorithm. 

With this scheme in place, only pseudo constraints will have multipliers that are out of bounds in the step 3 
optimality test. The number of non-optimal active constraints will never increase. In turn, each pseudo constraint 
will even tual ly be the constraint that gets chosen for dropping, though this may happen while another pseudo 
constraint is being relaxed. 

A special order has been chosen for relaxing pseudo constraints to take advantage of the possibility of savings 
from "stability". The modified main algorithm starts by testing and relaxing only stage-N pseudo constraints in steps 
3-6. This continues until all of the stage-N pseudo constraints have dropped from the active list or have zero 
multipliers. Then the algorithm switches to exclusive consideration of the stage-(N-l) pseudo constraints in steps 3-6. 
It performs add/drop cycles until all of these pseudo constraints get dropped or have zero multipliers. It continues this 
stagewise pseudo constraint relaxation scheme in a backwards sweep all the way to stage 0. The guessed solution is 
optimal after the last stage-0 pseudo constraint has been dropped or has had its multiplier go to zero. The trajectory 
from stage k to stage N is an optimal trajectory once all of the stage-k pseudo constraints have been dropped or have 
had their multipliers go to zero (although it probably will not be the final optimal trajectory associated with the 
solution to the overall problem). 

Comparison of Algorithm Complexity with Matrix Riccati Equation 

Table 1 compares the present algorithm's computational complexity with that of related algorithms for a typical 
aerospace controls problem. The time-varying multi-stage Matrix Riccati equation actually does not compute an A, ct 

because it solves a different optimization problem. It has been included because control engineers are more familiar 
with it. The three QL factorization entries assume that the factors are built up from initial pseudo constraints via 900 
rank-1 updates. In the last two entries, assumptions are made about the average number of stages affected per rank-1 
update. The table clearly indicates that the staircase QL factorization makes a tremendous improvement in comparison 
to the dense factorization; the improvement will not be nearly so great in comparison to a general sparse matrix code. 
Also, large improvements are expected from the special ordering of the pseudo constraint relaxation. Note that all of 
the algorithms are far more costly than the implementation of a time-varying LQR solution. Inequality constraints are 
difficult to handle. 

Table 1. 

A Comparison of Effort for Factorization of A ict for a Typical Aerospace Control Example 
(100 stages, 6 state vector elements, 3 control vector elements) 


Solution Method EflgH 

(No. of Mult., Div., & Sqrt) 


Matrix Riccati Equation 1 12,000 

Dense Matrix QL, Not storing Q 2,920,000,000 

Staircase QL, Arbitrary order of pseudo constraint relaxation 90,700,000 

Staircase QL, Special order of pseudo constraint relaxation 4,400,000 


Aerospace Example 


A simple aerospace control problem has been solved with the algorithm in order to demonstrate the usefulness 
of this class of algorithms on aerospace problems and in order to study the algorithm’s behavior. The problem is one 
of fixed-time maximization of the distance travelled by a thrust- and impulse-limited ground-to-ground rocket. The 
continuous-time problem is: 


find: u(t) for 0 < t < tf = 12 sec 

to minimize: J = - [ 1 0 0 0 ]x(tf) 

subject to: x(0) = 0 


(28a) 

(28b) 

(28c) 
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u + 


(28d) 


"0 1 00 -] 0 0 

dx_ 00 00 32.2 0 

dt - 0001 X+ 0 0 

-OOOoJ 0 32.2 


Hu(t)ll 2 < 5 g 



2 dx < lOg-sec 



(28e) 

(28f) 


-[OOlO]x(t) < Oft. (28h) 

which is a point-mass model of motion in the vertical plane. The acceleration (thrust) limit is 5 gs and the impulse 
limit is 10 g-sec. The first two state vector elements are horizontal position and velocity; the last two elements are 
vertical position and velocity. Only thrust and a uniform gravity field act on the rocket The first control vector 
element is horizontal acceleration; the second element is vertical acceleration. Constraint 28h keeps the rocket above 
the ground. 

In order solve this problem with this paper’s algorithm, the control time history has been approximated by a 
24-stage zero-order hold. Additionally, the norms in constraints 28e and 28f have been approximated by functions 
with octagonally-shaped contours. The 2-norm's contours are spherical; so, this approximation introduces some 
modelling error. Fixing the end time seems unnatural, but it is necessary in order to be able to model the problem as 
an LP. An NP model is needed to handle the fnee-end-time case. 

The LP code solved this problem in 55 min. on an IBM PC-AT with an 80287 coprocessor. It started from a 
first guess that violated inequality constraint 28h at every stage and that foolishly tried to maintain a constant thrust 
for the entire trajectory. Figures 1-3 compare the multi-stage LP solution with the exact continuous-time solution. 
In Fig. 1 , the LP solution does better than the exact solution because of mis-modeling; it takes advantage of some 
extra thrust available at some points of the octagon norm. The thrust magnitude and angle time histories. Fig. 2 & 3, 
are both close to the exact solution, and the discrepancies are due to the same modeling error. 

Figure 4 gives a 2-dimensional histogram of the constraint addition/deletion frequency. The left-hand horizontal 
axis indicates the stage at which the pseudo constraints are being relaxed in the special backwards-chaining process. 
The right-hand horizontal axis indicates the stage at which constraint additions and deletions are occurring during that 
relaxation process. The vertical axis gives the frequency of additions/deletions at the given right-hand-axis stage during 
pseudo-constraint relaxation at the given left-hand-axis stage. The extreme left-hand side of the figure shows no 
constraint addition or deletions — none can occur at any stage before stage k-1 when the pseudo constraints at stage k 
are the ones being relaxed. The peaks on and near the center diagonal of the graph lend support to the conjecture that 
most of the constraint additions/deletions will happen at stages near the pseudo-constraint-relaxation stage. Note, 
however, that a moderate amount of constraint addition/deletion activity occurred near the terminal stage throughout 
the optimization. Nevertheless, the average factor update was relatively cheap. Altogether, about 800 rank-1 updates 
occur during the optimization. The total number of decision vectors in the time history is 312 - 9 extra states are 
needed to model the impulse constraint in eq. 28f. 

Conclusions 

An algorithm has been presented for solving Dynamic Linear Programs. It takes advantage of the staircase 
structure of the active constraint matrix by factorizing it into staircase QL factors. These are derived in a stagewise 
fashion and play a role similar to that played by the time-varying matrix Riccati equation in multi-stage LQR theory. 
All of the usual linear programming functions have been implemented with the staircase QL factorization: decision 
vector solution, multiplier solution, and rank-1 updating. Each function has a computational complexity of 0(n 2 N) or 
less, where n is a block dimension and N is the number of stages. This is the best that can be expected for dense 
blocks. Numerical stability is assured via the exclusive use of QL factors and is independent of pivoting strategies. 

The algorithm is a modified active set implementation of the big M method with pseudo-constraint 
initialization. The modification restricts the set of non-optimal constraints that can be relaxed at one time to a single 
stage. This restriction gets iterated through all the stages in a backwards chain. Also, the modification chooses the 


339 




constraint that gets dropped in a way that assures optimality of the final portion of the solution time history. The 
modified strategy's goal is to reduce the average complexity of the rank-1 updates. 

A 24-stage example problem has been solved. The algorithm solves the 312-dimensional problem in about 800 
add/drop cycles, requiring 55 min. on an IBM PC-AT. The average update complexity is significantly reduced by the 
modified active set strategy. 
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Fig. 2 Thrust (Acceleration) Magnitude Time History for Ground- 
to-G round Missile Range Maximization 


Fig. 3 Thrust Angle Time History for Ground -to-G round Missile 
Range Maximization 



Fig, 4 Two-Dimensional Histogram of Frequency of Constraint 
Additions and Deletions at a Given Stage for a Given Pseudo- 
Constraint Relaxation Stage 





